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Most control engineering problems are characterized by several objectives, 
which have to be satisfied simultaneously. Two widely used methods for 
finding the optimal solution to such problems are aggregating to a single 
criterion, and using Pareto-optimal solutions. This paper proposed a Pareto- 
based Surrogate Modeling Algorithm (PSMA) approach using a combination 
of Surrogate Modeling (SM) optimization and Pareto-optimal solution to find 
a fixed-gain, discrete-time Proportional Integral Derivative (PID) controller 
for a Multi Input Multi Output (MIMO) Forced Circulation Evaporator 
(FCE) process plant. Experimental results show that a multi-objective, 
PSMA search was able to give a good approximation to the optimum 
controller parameters in this case. The Non-dominated Sorting Genetic 
Algorithm II (NSGA-II) method was also used to optimize the controller 
parameters and as comparison with PSMA. 
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1. INTRODUCTION 

In real control engineering world, engineers are often faced to track several objectives 
simultaneously. Most controllers are needed that provide a fast response, small overshoot, no oscillation and 
economical control. There are mainly two ways of tackling this problem: aggregating the objectives to a 
single objective or solving a multi-objective optimization problem using Pareto-based method. Aggregating 
several objectives into a single objective has the advantage of solving a simpler problem, but on the other 
hand many design iterations are required to obtain an acceptable compromise. On the other hand, the multi- 
objective approach is claimed to lead to a set of solutions each of which dominates the others in some sense. 

There are several of methods published and widely used to do multi-objective optimization for 
engineering problem such as NSGA and SPEA that are based on genetic algorithm and evolutionary 
algorithm. However despite of ability to achieve good optimization results [1], [2], both methods are known 
to need many function evaluations. In real engineering problem the cost of evaluating design is probably the 
biggest obstacle that prevents extensive use of optimization procedures. In the multi-objective world, this 
cost is multiplied, because there are multiple results to obtain. Evaluating directly a finite element model can 
take several days, which makes is very expensive to try hundreds or thousands design. Thus Pareto-based 
surrogate modeling algorithm (PSMA) is proposed for the determination of simpler models that involves less 
computational and gives good approximation results of the complicated model. 
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Surrogate Modeling (SM) also known as metamodeling or model reduction is said to be a model of 
a model or an approximation of a model. It is a supplementary model that can be alternatively used to 
interpret a more detailed model [3]. SM are usually consists of mathematical functions. These are functions 
with calibrated parameters, which are used as abstractions and simplifications of the simulation model [4]. In 
computer simulation, a SM is used to substitute a computationally expensive simulation model with a more 
efficient one. The basic idea of SM is to construct an approximate model using function values at some 
sampling points, which are typically determined using experimental design methods [5]. A SM exposes the 
system’s input-output relationship through a simple mathematical function [3]. Thus the simulation time for 
SM is less than that of the actual simulation model. 

Recently, as studied in [6], SM had been used to optimize various type of system, included the 
nonlinear system. Some of the systems that were successfully optimized using the SM technique are the 
Cartesian Coordinates Control of Hovercraft System [7] and the unmanned underwater vehicle [8], [9]. 
Through their study, they also had proved that the SM technique can optimize various types of controller 
parameters, for example, the fuzzy logic controller and the PID controller. 

The core of SM is a metamodel that gives the prediction of a system’s output. Although the output 
from metamodel is an approximate of actual measurement of complex model, it gives a good approximate of 
the actual value. The evaluation of output value is fast and provides enough information during design phase 
of a system [10]. Examples of metamodel are Radial Basis Functions Neural Networks (RBFNN), Kriging 
Models (KR), Polynomial Regression (PR), Multivariate Adaptive Regression Splines (MARS), and Support 
Vector Machines (SVM). In comparison, RBFNN shows a generally better performance. Based on different 
types of problems (i.e., different orders of nonlinearity and problem scales) it is concluded that RBFNN is the 
most dependable method in most situations in terms of accuracy and robustness [11]. In this project, a 
RBENN was used as the metamodel to approximate the mapping of the controller gains and the objective 
function. 


2. MODELING OF THE SYSTEMS 
2.1. Radial Basis Function Neural Network 

Radial Basis Function Neural Network (RBFNN) was used as the Metamodel to approximate the 
mapping of the controller parameters and the objective function. The radial basis functions were first used to 
design Artificial Neural Networks in 1988 by Broomhead and Lowe [12]. The architecture of the RBF NN 
used in this work is illustrated in Figure 1. 


Hidden layer 


Input layer 


Figure 1. Radial Basis Function Neural Network 


The network consists of three layers: an input layer, a hidden layer and an output layer. Here, R 
denotes the number of inputs while Q the number of outputs. Equation (1) is used to calculate the output of 
the RBF NN for Q = 1, the output of the RBFNN in Figure 1 is calculated according to 


n(x, w) = Lud =C I), 


(1) 


Rx1 
Where ¥ © R is an input vector, ( Jig a basis function, | l denotes the Euclidean norm, 1% 
are the weights in the output layer, S1 is the number of neurons (and centers) in the hidden layer and 


Rxl 
eR ; ; ; i 
Ck are the RBF centers in the input vector space. Equation (1) can also be written as Equation (2) 
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19(x,w) =" (x)w 


(2) 
Where basis function in Equation (3) 
P (x)=[4(|z—al)---4s (e-esil) 6) 
And weight layer in Equation (4) 
w” =[w; mz- -ms] (4) 


The output of the neuron in a hidden layer is a nonlinear function of the distance given by Equation (5): 


_ ee 
(x)=e (5) 


Where B is the spread parameter of the RBF. For training, the least squares formula was used to find 
the second layer weights while the centers are set using the available data samples. Thus, the approach of 
Pareto-based Surrogate Modeling Algorithm (PSMA) for multiobjective optimization as summarized in [6-9] 
was used in this project. 


2.2. Forced Circulation Evaporator 

In addition, a metamodeling approach for PID controller in an evaporator process has been 
successfully presented in [13], [14]. Figure 2 shows the forced circulation evaporator derived by Newell and 
Lee [15] in 1989. This evaporator has become a well-known and very difficult benchmark used by control 
engineers to evaluate their methodologies. A feed stream enters the evaporator with concentration X1, 
temperature T1 and flow rate F1. It will mix with recirculation liquor, which is pumped through the 
evaporator at flow rate F3. The evaporator itself is a heat exchanger, which is heated by steam flowing at a 
rate F100, with temperature T100 and pressure P100. The mixture of feed and recirculation liquor boils 
inside the heat exchanger, and the resulting mixture of vapor and liquid enters the separator, which the liquid 
level is L2. The operating pressure inside the evaporator is P2. Some portion of liquid from separator drawn 
out as product with concentrationX2, with flow rate F2 and temperature T2; most of it becomes the 
recirculation liquor with flow rate F3. The vapor from the separator flow to a condenser at flow rate F4 and 
temperature T3, where it is condensed by cooled water flowing at flow rate F200, with entry temperature 
T200 and exit temperature T201. 
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Figure 2. Forced Circulation Evaporator 
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The constant value and description are shown in Table 1, while the variables names, descriptions, 
steady state value, and engineering units are shown in Table 2. 


Table 1. Constant value and description 


Constant Description Value Units 
pA Liquid density and cross-sectional area of separator 20 kg/m 
M Amount of liquid in the evaporator 20 kg 

C Constant that converts the mass of vapor into an equivalent pressure 4 kg/kPa 
Cp Heat capacity of the liquor 0.07  kg/min 
A Latent heat of vaporization of the liquor 38.5 kg/min 
As Latent heat of steam at the saturated conditions 36.6 kg/min 
Uaz Overall heat transfer coefficient times the heat transfer area 6.84 kW/K 


Table 2. Evaporator variables and steady state value 


Variable Description Value Units 
F; feed flow rate 10.0 kg/min 
F, product flow rate 2.0 kg/min 
F circulation flow rate 50.0 kg/min 
Fy vapor flow rate 8.0 kg/min 
F; condensate flow rate 8.0 kg/min 
Xi feed composition 5.0 percent 
Xo product composition 25.0 percent 
Tı feed temperature 40.0 deg C 
To product temperature 84.6 deg C 
T; vapor temperature 80.6 deg C 
L separator level 1.0 metres 
P2 operating pressure 50.5 kPa 

Fioo steam flow rate 9.3 kg/min 
Tioo steam temperature 119.9 deg C 
Pioo steam pressure 194.7 kPa 
Qioo heater duty 339.0 kW 
F00 cooling water flow rate 208.0 kg/min 
T200 cooling water inlet temperature 25.0 deg C 
Toor cooling water outlet temperature 46.1 deg C 
Q200 condenser duty 307.9 kW 


3. RESEARCH METHOD 

Two control variables are chosen out from FCE as objectives function and controlled by using PID 
controller. The control variables are L2 and P2 with manipulated variables of the plant are F2 and F200. The 
design objective will be a six parameter optimization problem of determining the optimal parameter gains 
[Kp1 Kil Kd1 Kp2 Ki2 Kd2] to minimize the output of L2 and P2. Table 3 shows parameter coefficient with 
their range which cover both PID controllers. This range is used in order to obtain a good comparison 
between PSMA and NSGA-II. 

The simulation for both controllers was done using MATLAB SimulinkTM as illustrated in 
Figure 3. All values were initialized at the operating points as stated in Table 2. Simulation time was set to be 
300 seconds and run using ode14x (extrapolation) solver. The set point for P2 is 50.5 kPa over the simulation 
time while L2 was given a varying step input from initial 1.0m to 2.5m and going down back to 1.0m. 
Table 4 shows the control variable constraints. 


Table 3 PID variables and design space 
Variables 
Limit Kyi K; Ka Ky2 K; Kaz 
Lower -130 -2 -60 -410 -20 -10 
Upper -100 2 -50 -390 -10 -5 


Table 4 Variable constraints 


Variable Lower limit Upper limit 
F 0 kg/min 50 kg/min 
F200 0 kg/min 400 kg/min 
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The performance criterion to measure the output tracking in this case was the Integral Square Error 
(ISE) given by: 


2 
ISE = [(y, (t)- y(t)?’ dt 5 
Where yg is the desired output (set point) while y is the actual output. This criterion has been used 


because of the ease of computing the integral both analytically and experimentally. The most efficient value 
of Pareto frontier is defined by calculating Euclidean distance between ISE and initial point, zero: 


Cost = [S ase? 
i=l 


Figure 3 shows PID Forced Circulation Evaporator as implemented in Matlab® Simulink®. 
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Figure 3. PID Forced Circulation Evaporator as implemented in Matlab® Simulink® 


3.1. Pareto based Surrogate Modeling Algorithm 


Table 5 show the objective function, initial design space (D) and larger design space (D’) used for 
PSMA simulation. 


Table 5. Objective function, initial data sets and large data sets 
PID 


Objective function Initial data sets (D) Large data sets (D’) 


Parameter 

Kyi {-130, -120, -110} {-130, -125,..., -110} 
F2 Kir {-2, 0, 2} {-2, -1,..., 2} 

Ka {-60, -55, -50} {-60, -55, -50} 

K,2 {-410, -400, -390} {-410, -405,..., -390} 
P2 Kz {-20, -15, -10} {-20, -17.5,..., -10} 

Ka {-10, -5} {-10, -7.5, 5} 

Total number of data configurations 486 5625 


The step size of D and D’ specifically sets by user where D’’ use smaller resolution thus multiplies 
the total number of data configuration. Different with NSGA-II, the value between bound are created 
randomly. The initial data sets should not too small for proper training and should not be too large to 
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minimize the training time. The initial data sets are used to simulate ISE for both operating pressure P2 and 
separator level, L2 simultaneously. RBFNN then use ISE value from initial data sets and predict the output 
for large data sets. 

In this PSMA simulation, the basis function centers, c, is set equal to the input vector from the 
training set or maximum number of initial data sets, 486. The spread value of 10 is used in the training 
process. The larger the spread of the data the smoother will be the function approximation. A large spread 
implies a lot of neuron will be required to fit a fast changing function. Where a small spread is means less 
neuron will be required to fit a smooth function and the network may not generalize well. 


3.2. Non dominated Sorting Geneting Algorithm 

The NSGA-II [16] is selected as comparison to PSMA because of widely used and capable 
algorithm. The principle behind NSGA-II is that the non dominated solution that usually occur for 
multiobjective optimization problems are all treated as equals. This allows the algorithm to evolve a set of 
non-dominated solution that is equally well suited for solving the specific problem given the performance 
measures specified. By using the algorithm for tuning of PID controller for the FCE, it will be possible to 
obtain varied set of different solution that should perform well with regards to minimization of all specific 
performance measures. NSGA-II run-time parameters used for this problem are summarized in Table 6. 

The choice of real valued representation was made to ensure that the precision of the parameters 
would not be compromised by a choice of precision, which can happen for binary representation. A crossover 
probability of 0.9 ensures a good mixing of genetic material and mutation probability can be expressed as 


1 


param where = ?4"“"" is the number of parameters in an individual which for this application is six. 


Simulated binary crossover parameter (SBX) and the mutation parameter were decided to use 20 and 20 
respectively since they provide a reasonable distribution of solutions for the different operations. 


Table 6. NSGA-II run-time parameters 


Representation type Real values 
Crossover probability 0.6 
Mutation probability 0.167 
SBX parameter 20 
Mutation parameter 20 
Population 100 
Generation 100 


4. RESULTS AND ANALYS 
4.1. Simulation Result of PSMA 

Figure 4 show the simulation result of P2 and L2 using initial data sets with 486 total number of 
data configurations. 
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Figure 4. ISE of separator level and operating pressure for initial data sets. 
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The ISE values are used to train the RBFNN which will then be used as the metamodel of the FCE 
to evaluate the ISEs for the corresponding large data sets of the controller parameters. The results of RBFNN 
training using 486 centers and 10 spread are shown in Figure 5. 

After the training stage RBFNN is used to perform the simulation for large data space controller 
parameters sets which consist of 5625 data sets. The result is shown in Figure 6. The estimated ISE for L2 
and P2 then plotted into pareto set as in Figure 10. The pareto-optimal frontier marked with blue circle. Since 
the both objective function to find minimum value, the closest to origin indicates the most efficient value, 
represented by green triangle in the figure. Although the most efficient value predicted by RFBNN (5.417 for 
L2 and 6.744 for P2) is not same with real simulation, PSMA was able to give minimum coefficient 


parameter as in Table 5. 
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Figure 5. RBFNN training of ISE of initial data sets 
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Figure 6. Surrogate modeling output for large data sets of L2 and P2 
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Figure 7. Plot of Pareto optimal frontier for L2 and P2 using PSMA 
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4.1. Simulation Result of NSGA-II 

As comparison to PSMA, Figure 8 show the Pareto set of NSGA-II optimization. Most efficient 
value marked with green triangle. 

Parameters of PID controller and their relevant cost values obtained by PSMA and NSGA-II 
approach are demonstrated in Table 7. From the simulation results in Table 8, the parameter controller 
obtained by PSMA clearly has better performance than NSGA-II. The ISE value obtained by PSMA for both 
outputs, L2 and P2 is lower than using NSGA-I. The PSMA simulation time took only 1.52 minutes 
compare to NSGA-II, 23.36 minutes. In general the controller obtained by PSMA has the best performance. 
The result in Table VIII shows the ability of PSMA in dealing with challenging optimization problems. 

Figure 9 shows the response of controlled FCE to step input using different controllers obtained by 
PSMA and NSGA-II. 


Multiobjective Optimization Using NSGA-II 
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Figure 8 Plot of Pareto optimal frontiers for L2 and P2 using NSGA-II 


Table 7. Parameter of PID controller obtained by Table 8. ISE, Cost and simulation time by PSMA 


PSMA and NSGA-II and NSGA-II 
Method Method 

PID Parameters PSMA _NSGA-II Criteria PSMA NSGA-II 
Ka -125 -119.66 ISE for L2 8.041 8.187 
Ka 2 1.89 ISE for P2 6.618 6.801 
Kar -55 -59.05 Cost 10.410 10.640 
Kp2 a0 -405.61 Simulation time (min) 1.520 23.360 
Ko -20 -19.14 
Ko -10 -9.61 
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Figure 9 (a) Response of separator level, L2. (b) Response of operating pressure, P2 
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The controllers gave a good response for separator Level, L2. In Figure 9(a) the settling time by 
using PSMA parameter is slightly better than NSGA-II. It can be seen at second 250 when step input changed 
to set point 1m, PSMA respond reach steady state until second 300. For operating pressure, P2, response 
obtained by NSGA-II and PSMA parameter almost identical. This condition occurs because the parameter 
gains, Kp2 Ki2 Kd2 of both optimization technique almost the same. 

Similar to other optimization algorithm such as SPEA, NSGA, the discussed method in this paper, 
PSMA does not necessarily guarantee the real time requirements in exact applications. But as shown in this 
paper, PSMA was able to give fast computational time to obtain best value for the controller. In application 
of high computational complexity, the use of PSMA will be more preferable. 


5. CONCLUSION 

The purposed optimization method using PSMA offers advantages at especially reducing the cost 
and time by utilizing surrogate modeling for complex and expensive design. The genetic algorithm based 
optimization required large number of objective function evaluation to generate Pareto-optimal front. 
Therefore the evaluation of the required number of objective function values through a full model 
experiment. In this study NSGA-II took around 15 times simulation time to optimize the operating pressure 
and separator level of FCE whereas PSMA training and testing takes couple of minutes depending upon the 
user’s experience and prediction through surrogate modeling. The PSMA approach us clearly a useful 
approach and this will become more significant for a larger D of for a more complicated problem. 

Using FCE as a study case, PSMA used to optimize the parameter gain of PID controller. Surrogate 
modeling does provide the designer with a quick estimate for a good set of good parameter to begin with. 
Further simulation on the actual system can be done if better values are required. In this example, the data set 
D was created by choosing the input values like the grid fashion based on background knowledge of the 
problem. A more intuitive approach is to start with a small number of samples and then sequentially add 
more data samples intelligently employing Experimental Design techniques such as Worst Case Approach 
and Cross Validation technique. It is envisaged that a more strategic data location will allow the creation of a 
more accurate surrogate modeling using less data, therefore, less time required to estimate the best controller 
parameters. 
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